home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / ODEINT.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  81 lines
  1. PROCEDURE odeint(VAR ystart: glnarray; nvar: integer;
  2.        x1,x2,eps,h1,hmin: real; VAR nok,nbad: integer);
  3. (* Programs using routine ODEINT must provide a
  4. PROCEDURE derivs(x:real; y:glnarray; VAR dydx:glnarray);
  5. which returns the derivatives dydx at location x, given both x
  6. and the function values y. They must also define the type
  7. TYPE
  8.    glnarray = ARRAY [1..nvar] OF real;
  9. and must declare the following parameters
  10. VAR
  11.    kmax,kount: integer;
  12.    dxsav: real;
  13.    xp: ARRAY [1..200] OF real;
  14.    yp: ARRAY [1..10,1..200] OF real;
  15. and must initialize kmax and dxsav in the main routine. *)
  16. LABEL 99;
  17. CONST
  18.    maxstp=10000;
  19.    two=2.0;
  20.    zero=0.0;
  21.    tiny=1.0e-30;
  22. VAR
  23.    nstp,i: integer;
  24.    xsav,x,hnext,hdid,h: real;
  25.    yscal,y,dydx: glnarray;
  26. BEGIN
  27.    x := x1;
  28.    IF (x2 > x1) THEN h := abs(h1) ELSE h := -abs(h1);
  29.    nok := 0;
  30.    nbad := 0;
  31.    kount := 0;
  32.    FOR i := 1 TO nvar DO BEGIN
  33.       y[i] := ystart[i]
  34.    END;
  35.    IF (kmax > 0) THEN xsav := x-dxsav*two;
  36.    FOR nstp := 1 TO maxstp DO BEGIN
  37.       derivs(x,y,dydx);
  38.       FOR i := 1 TO nvar DO BEGIN
  39.          yscal[i] := abs(y[i])+abs(dydx[i]*h)+tiny
  40.       END;
  41.       IF (kmax > 0) THEN BEGIN
  42.          IF (abs(x-xsav) > abs(dxsav)) THEN BEGIN
  43.             IF (kount < kmax-1) THEN BEGIN
  44.                kount := kount+1;
  45.                xp[kount] := x;
  46.                FOR i := 1 TO nvar DO BEGIN
  47.                   yp[i,kount] := y[i]
  48.                END;
  49.                xsav := x
  50.             END
  51.          END
  52.       END;
  53.       IF (((x+h-x2)*(x+h-x1)) > zero) THEN h := x2-x;
  54.       rkqc(y,dydx,nvar,x,h,eps,yscal,hdid,hnext);
  55.       IF (hdid = h) THEN BEGIN
  56.          nok := nok+1
  57.       END ELSE BEGIN
  58.          nbad := nbad+1
  59.       END;
  60.       IF (((x-x2)*(x2-x1)) >= zero) THEN BEGIN
  61.          FOR i := 1 TO nvar DO BEGIN
  62.             ystart[i] := y[i]
  63.          END;
  64.          IF (kmax <> 0) THEN BEGIN
  65.             kount := kount+1;
  66.             xp[kount] := x;
  67.             FOR i := 1 TO nvar DO BEGIN
  68.                yp[i,kount] := y[i]
  69.             END
  70.          END;
  71.          GOTO 99
  72.       END;
  73.       IF (abs(hnext) < hmin) THEN BEGIN
  74.          writeln('pause in routine ODEINT');
  75.          writeln('stepsize too small'); readln
  76.       END;
  77.       h := hnext;
  78.    END;
  79.    writeln('pause in routine ODEINT - too many steps'); readln;
  80. 99:   END;
  81.